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Abstract 

We investigate asymptotic decay phenomenon towards the nonequilibrium steady 
state of the thermal diffusion in the presence of a tilted periodic potential. The pa- 
rameter dependence of the decay rate is revealed by investigating the Fokker-Planck 
(FP) equation in the low temperature case under the spatially periodic boundary 
condition (PBC). We apply the WKB method to the associated Schrddinger equa- 
tion. While eigenvalues of the non-Hermitian FP operator are complex in general, 
in a small tilting case the imaginary parts of the eigenvalues are almost vanishing. 
Then the Schrddinger equation is solved with PBC. The decay rate is analyzed in 
the context of quantum tunneling through a triple-well effective periodic potential. 
In a large tilting case, the imaginary parts of the eigenvalues of FP operator are 
crucial. We apply the complex-valued WKB method to the Schrddinger equation 
with the absorbing boundary condition, finding that the decay rate saturates and 
depends only on the temperature, the potential periodicity and the viscous con- 
stant. The intermediate tilting case is also explored. The analytic results agree well 
with the numerical data for a wide range of tilting. 
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1 Introduction 



Thermal diffusion process in a periodic potential in the presence of external 
force appears in many situations. In solid-state physics, it appears in diffu- 
sion of ions or molecules on crystal surfaces [1], Josephson junctions [2], motion 
of fluxons in superconductors [3], rotation of molecules in solids [4], superionic 
conductors [5], charge density waves [6], to mention a few. As a ratchet system, 
it is also applied for particle selection[7], biophysical processes[8] and intra- 
cellular transport [9] [10]. Most of the researches in this context are concerned 
with the steady state. For example, it is shown that a symmetric dichotomous 
noise causes non-zero steady state current in such a system[ll]. The existence 
of an optimal diffusion coefficient vs. external force is also reported [12]. 

Such a diffusion process can be described by the over-damped Langevin equa- 
tion with a tilted periodic potential, 



<*(*)> = o, x mm) = w5(t-t'). a) 

Here, r], — ksT and W are the viscous coefficient, temperature and the 
external force, respectively. is the Langevin force which satisfies the fluc- 
tuation dissipation theorem. The periodic potential U^\x + L) — U^\x) 
with period L is tilted by the external tilting force F = ^j^- and we call the 
potential U (x) = (x) — '^-x as a tilted potential. 

In order to investigate the probability density P(x, t) that the particle is found 
in a position x at time t, it is useful to transform the Langevin equation into 
the corresponding Fokker-Planck(FP) equation, 



Hereafter we solve FP equation (2) under the (spatially) periodic boundary 
condition (PBC), P(x+L, t) = P(x, t), which is satisfied by the steady state at 
t — > oo. The analytic expression for the steady state P st (x) = lim^oo P(x, t) 
under PBC is well-known[9,13,14]: 



P°\x) = Ny U -^ X fdye^ ) 

X 
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(3) 



where N is the normalization constant chosen so that Jq dxP st 
steady state has a non-zero probability current J st (x) given by 



(x) = 1. The 



J*t = l( U(x) - 9-pi P s \x) — N (l — 



v V 
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The above steady state is a nonequilibrium state. In this paper, we study the 
low-temperature relaxation dynamics towards the steady state in the presence 
of a tilted periodic potential. In particular, asymptotic behaviors of the decay 
phenomenon is investigated by the eigenfunction expansion method. Our main 
interest lies in parameter dependence of the decay rate. 

This paper is organized as follows. In section 2 we present the framework 
of our analysis, relevant numerical data for eigenvalues of the FP operator, 
and the main result. Section 3 treats the important theme on the boundary 
condition for eigenstates of the associated Schrodinger equation. In subsection 
3.1 the small tilting case is investigated in the context of a tunneling problem. 
In subsection 3.2 the large tilting case is treated, where we shall apply the 
complex-valued WKB method to the Schrodinger equation with the absorbing 
boundary condition deduced from PBC for the original FP equation. Section 
4 is devoted to a summary. Appendix is concerned with detailed calculations 
in section 3.2. 



2 Eigenvalues of the Fokker-Planck operator 

The WKB treatment of FP equation has a long history since the monumental 
work by van Kampen[15]. Most of the studies, however, were limited to the 
system with a well-defined thermal equilibrium state[16]. This fact guaran- 
teed a real-valued nature of eigenvalues for the non-Hermitian FP operator, 
and one can reduce FP equation to the eigenvalue problem of the associated 
Schrodinger equation with the same natural boundary condition as the FP 
equation has. On the other hand, relaxation dynamics towards the nonequi- 
librium steady state with the non-vanishing constant current demands a qual- 
itatively different WKB approach. In fact FP operator in such a system has 
complex eigenvalues in general. 

This situation is recognized in our system with a tilted sinusoidal potential [17], 




(5) 
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where Uq is the amplitude of the periodic part of the potential. Putting eigen- 
values of the non- Hermit ian FP operator as — E and substituting P(x, t) = 
e~ Et P(x) into Eq.(2), we have the eigenvalue problem, 



While the steady state in Eq.(3) is periodic, P st (x + L) = P st (x), and is char- 
acterized by the zero Bloch wavenumber {k — 0), it is difficult to analyze the 
relaxation dynamics in general. Any initial distribution P(x, 0) which include 
k 7^ Fourier components can relax towards P st (x): the fc ^ components 
will decay out on the way of time evolution. In the present work, however, we 
confine ourselves to the manifold of k = 0, which makes the analysis of the 
decay phenomenon accessible. This implies that the time-dependent solution 
of FP equation in Eq. (2) relaxing towards the steady state should satisfy the 
spatial PBC, P(x + L,t) = P(x, t) and be constructed from solutions with a 
natural boundary condition Q(x,t) (lmx E _ > ± 0O Q(x, t) = 0) as 

oo 

P(x,t)= Yl Q(x + nL,t). (7) 

n=— oo 



Since the FP operator in Eq.(2) is invariant against space translation by L, it 
has Bloch-type eigenstates {P nt k{x)} and (complex) eigenvalues {E n ^}- Noting 
the normalization matrix J L dxP£ n (x)Pk' , n '(x) = ^k,k'N^ n , due to the non- 
Hermitian nature of FP operator, the solution P(x, t) in Eq.(7) is now written 
explicitly as 

P(x, t) = J2 C n ,k=oPn,k=o(x)e~ En ^°\ (8) 

n 

L 

C n , k=0 = E(^ (0) )ni / dxP:,, k=0 (x)P(x, 0), Cn^o = 0, (9) 

n' o 

which starts from an arbitrary k = distribution P(x, 0) at t — and relaxes 
towards P st (x). 

In the manifold with k — 0, the probability distribution can be expanded as 
P{x) = J2n c n,e^ rn - Then the eigenvalue problem (6) is reduced to 



nUn 



(10) 
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Fig. 1. Distribution of the complex eigenvalue with the second smallest real part 
for full parameter range 3 < Uq < 20, < W < 20 (dots), and for the case of fixed 
potential barrier Uq = 10, < W < 20 (solid line). A pair of solid lines rapidly 
merge to the real axis as W decreases below Uq. 

which can be solved numerically by truncating the infinite dimensional matrix 
A to a large but finite dimendional matrix. 

The zero-eigenvalue (E = 0) corresponds to the steady state. Under PBC, as 
is obvious from Eq. (8), the probability distribution approaches exponentially 
to the unique steady state, and the decay rate A is obtained as the real part of 
the second smallest eigenvalue. Among many complex eigenvalues of Eq.(10), 
therefore, we have focused on the one with the second- smallest real part and 
its variation against various parameter values Uq,W (see Fig.l). From Fig.l 
one can see that the eigenvalues with large positive real parts have nearly zero 
imaginary part. (See tails of solid lines in Fig. 1.) These eigenvalues correspond 
to the small tilting case. In fact, a pair of solid lines rapidly merge to the real 
axis as W decreases below Uo, ensuring the real-valued nature of eigenvalues 
in the case of small tilting. 

Uo and W dependence of the decay rate A is shown in Fig. 2, which we have 
confirmed by numerically solving the FP equation in Eq.(2). There is a smooth 
but steep decrease of the decay rate around the crossover line of Uo ~ W, 
where the original potential minima disappear. 

In the next section, the Schrodinger equation associated with FP equation 
will be solved with use of the WKB method. In the small tilting case in 
section 3.1, we impose PBC for the wavefunction of the Schrodinger equa- 
tion and apply the Bloch theorem. The decay rate for W < Uo is shown to 

be [^) \Jl — (^) Uq- In the large tilting case in subsection 3.2, we have 
recourse to the absorbing boundary condition and apply the WKB method 
which admits complex energies. The decay rate for W ^> Uo is found to be 
(2it/L) 2 6 (with no W and U dependence). The crossover region where W is 
comparable to Uq is numerically analyzed based on this approach. 
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Fig. 2. W and Uq dependence of the decay rate A. 
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Fig. 3. W dependence of the decay rate A for Uq = 3, 9 = 0.1. Dots show the decay 
rates obtained from numerical diagonalization of FP operator. Three solid lines are 
theoretical results (For small tilting case W < Uq the decay rate A is obtained as 
the first excited level huj a of the associated Schrodinger equation (left part). For 
crossover regime W ~ Uq, A is given as the numerical solution of (37) (middle part). 
For large tilting case W S> Uo, the perturbative solution of (38) gives A (right part). 

3 WKB analysis of probability distribution 

The FP equation is transformed into the Schrodinger equation by the separa- 
tion ansatz[15] 

P(x,t) = e-^r (j)(x)e- Et . (11) 
Then one has 

9^-<i){x) + {E-V{x))cf>{x) = U, (12) 
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where we redefined eigenvalue as Er] — > E, and defined the effective potential 



v(x)se{{ m)'.m } . (13) 



The Schrodinger operator (SO) in Eq.(12) looks Hermitian. Noting in Eq.(f f) 
U(x+L) = U(x)—2irW and P n ^(x+L) = e tkL P n ^(x), however, the wave func- 
tion 0(x) of the Schrodinger equation should satisfy the absorbing boundary 
condition, 



<p(x + L) =e ikL ~^(f)(x), (14) 



2ttW 



which eventually breaks the Hermitian nature of SO. 

In the following we shall analyze the Schrodinger equation in Eq.(12) first in 
the small tilting case and then in the strong tilting case. In the small tilting 
case, the numerical evidence in the previous section has shown the real- valued 
nature of the second-lowest eigenvalue, indicating the recovery of Hermitian 
nature of SO. Therefore we reduce Eq. (14) to 

(j)(x + L) = e lkL (j)(x), (15) 



by assuming the absorbing term being vanishing. The validity of this 
assumption will be verified by comparing between the analytic and numerical 
decay rates. Below we shall confine to the k = manifold. 



3.1 Small tilting case 



Substituting the tilted sinusoidal potential (5) into (13), one has the effective 
potential, 

vr 2 / . 2nx „ A 2 1 /2tt\ 2 „ 2tt 



V( X ) = m ^„sin— + W) + - [ T ) U cos TX . (16) 



Since the first term is much larger than the second one in the low temperature 
regime 9 U , the second one is omitted hereafter. 

In solving Eq.(12) the low temperature condition allows us to use the WKB ap- 
proximation, since 9 corresponds to ^ in the standard Schrodinger equation. 
In the case that the external force is weak (W < Uq), the effective potential 
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Fig. 4. Effective potentiai V(x) and the real part of the eigenvalue of FP operator 
with second smallest real part E with the corresponding classical turning points 
Xq, x 5 for W = 1, Uq = 3. 
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Fig. 5. Original tilted potential f/(x) in real space. Small tilting case (a), crossover 
region (b), and large tilting case (c). 

V(x) has three wells. We consider the energy levels for the eigenstates which 
experience all the three wells[18] (see Fig. 4). 

In terms of the tilted potential U(x), the weak external force guarantees the 
existence of potential minima in Fig. 5. Let us denote the wave functions in 
each region a, • • -, e bordered by classical turning points xq, • • •, x$ as <f> a , • • <p e 
and introduce the action integrals 



S a = - J pdx, M b = - J pdx 

Sc= hl pdx ' Md= hl pdx ' (17) 



X2 X 3 



where h 2 = 2m6 and p = ^2m\E — V(x)\ (p stands for the momentum in a 
potential- well region). S and M are defined in each well and barrier, respec- 
tively. Each neighboring wave functions are mutually related by the connection 
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formula 



J_ e ±i(s+5) (£ > 1 (> ± i e -s ) ( £ < l /( a; )) ) ( 18 ) 

VP \'l> v - - 



where 



2/0 

S=\-Jpdx\, (19) 



with |/o being the classical turning point. 
For example, <p a and 0^ are related as 



i i \ 1 / — t r X P<&; , t pdx 

> a {x) = — Cie ftJ *o +c 2 e hJx o 
VP V 



= JL Jde"* (e*£ , " fa + i e ^A^) (20) 

+C2e <s ^ /; **** _ y-k /; ^ j (21) 

In the same way, one obtains the expression for <f> e (x) as 

e (x) = |(A Cl + Bc^e^^ + (C Cl + Z}c 2 )e^>*} , (22) 

where A,B,C and .D are given as 

A = e iS ° {i (ie M » sin 5 C + cos 5 C ) e - Md + 2 (e Mb cos S c + ^e"^ sin s)j e Md | , 

5 = e - iS " |i (e M > sin S c + ^e- Mb cos 5 C ) e~ Md - 2 (ie^ cos S c + sin 5 C ) e Md | , 

C = e iS ° | X - [e Mb sin S c - l -e~ Mb cos 5 C ) e~ Md - 2 (-ie M * cos S c + \e~ Mb sin 5 C ) e Md } 

D = e~ iS " |i (-*e Mi > sin S c + ^e"^ cos 5 C ) e- Md + 2 (e M > cos S c - V M » sin S c ) e Md 

(9? 



(23) 



On the other hand, with use of the Bloch theorem in Eq.(15) based on the 
periodicity of the effective potential V, (f> e and (f> a should be related as 



e = e* fcL a , (24) 
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where k is the wave number. This relation gives a constraint to guarantee the 
nontrivial coefficients (ci, c 2 ) of (p a and e : 



A - e ifci B 
C D- e ikL 



0, 



(25) 



which can be rewritten as 



e ^ = COs{k L )+ iM k L ) = 1 -(A + D ) ± i J 1 - 1 - { A + D y. (26) 



Here we have used the fact AD — BC — 1. At low temperatures the tunneling 
integrals are sufficiently large, Mb ^> 1 and Md 1. Then we can make the 
approximation: 



■(A + Df = 4e 2M "e 2Md cos 2 (,S a ) cos 2 (S c 



(27) 



Since (27) should satisfy (26), | cos 5^1 and | cos S c \ must be sufficiently small 
and be approximated as 



COs(S a , c ) ~ (-1)"- {S a , c - K, c + 1)} 7T 

sin(S 0)C )~(-l) n <- (28) 

Furthermore, using the harmonic approximation in the neighborhood of the 
potential minima x^ in , the actions can be evaluated as 

S a ,c = ^-{E-V 0a>c ), (29) 



where tu a>c = \j j^V" (x a ^ in ) is the curvature of the potential minima, and V 0ajC 
stands for the local minima of V(x). Then Eq.(26) gives the quantization 
condition for eigenvalues E 



cos(fcL) = -(A + D) 



- | (^ e Mb+Md + i e " Mi,_Md 
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Solving the quadratic equation of E in (30), we obtain the lowest few eigen- 
values within the k — manifold, 



E± — E 0+ ± yj Eq_ + e 
~E 0+ ± \E - + 



(31) 



where E ± stand for the harmonic oscillator levels at the wells and e is a 
correction due to tunneling, which are defined as 



E ± = ^ <jV 0a ± Vqc + hw a [n a + ^ ± hj c (n c + ^ J , 

, h 2 u a uj c e~ 2M > + e~ 2M * + 2(-l) n « +n c e - Mb - M * cos /c 
e = e (k = 0) = — 1 _ 9A/f _ 9Mj . (32) 



71-2 4 _|_ lg-2M 6 -2JW d 

In general, the eigenvalues form subbands with their width given by the tun- 
neling term e, although we confine ourselves to the k — state in each of 
the subbands. At x a min = -^(niv — (— l) n arcsin ^), we find V 0a = V(x min ) ~ 



1 



(33) 



Noting further that Vo a + \hw a ~ Vo c + — 0, one can confirm the lowest 
eigenvalue E_ ~ 0, which is responsible to the steady state. Since the tunneling 
term is small compared to hu a , one finds the decay rate (:the second- lowest 
eigenvalue), 



A = E + ~ hhj a 



1 - 



w \u 



(34) 



We have confirmed numerically that this expression gives a good approxima- 
tion for the decay rate in the small tilting case (see Fig. 3). 

Now we can explain the W (external force) dependence of the decay rate. 
From (34) one has the monotonic W (external force) dependence of the decay 
rate. Intuitively, this parameter dependence of the decay rate is explained in 
terms of tunneling process as follows. The energy level ftw a = ^IQV'ix^^ 
is a monotonic function of the curvature of the potential minimum V"{x a mi ^) 
which decreases for increasing W, because the potential barrier becomes lower 
and finally the barrier disappears as W increases. The case that each barrier 
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separating adjacent local potential minima vanishes is treated in the next 
section. 



3.2 Large tilting case 



The separation ansatz P(x,t) = e w <p(x)e~ Et in Eq.(ll) seems to contain 
a divergent term coming from e^^sr involved in e ^~ . This problem was 
overcome by choosing the absorbing boundary condition in Eq.(14). Confining 
ourselves to k = states as in the previous subsection, we have 



4>{x + L) = e 29 (j)(x). (35) 

With use of the complex-valued WKB method, Eq. (35) leads to the quanti- 
zation condition for the complex eigenvalues E = E 1 + iE 2 as 

with V(x) given in Eq.(16). 

Remarkably there is no classical turning point on the real axis, and the energy 
becomes complex, as is confirmed numerically for the strongly tilted case (see 
Fig.l). Then the decay rate A is given by the real part of the eigenvalue with 
the second smallest part, which we denote E\. The existence of the imaginary 
part, which we denote E 2 , implies that the asymptotic behavior is oscillating. 

Noting again that the second term of (16) is suppressed at low temperatures, 
the WKB quantization condition (36) is rewritten as 



/ \ 1 + 2— sin + — s n 2 — -(E 1 + iE 2 ) 

9L J V W L \WJ L ttW 2V ; 



= — + %2ixn. (37) 

V 

For the large tilting case W 3> U , we can solve this equation perturbatively 
with the expansion parameter U /W (see Appendix). Since it is numerically 
shown that the integer n is ±1 for the parameter range 3 < Uq < 20, Uq < 
W < 20, we investigate the case that n = ±1. As a result, we obtain for large 
W 
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(38) 



These results are consistent with those from the numerical diagonalization of 
FP operator. It is interesting that for large W , the decay rate E\ saturates and 
depends only on the temperature, the period and the viscosity. It is also worth 
noting that E 2 , which represents the oscillation frequency in the asymptotic 
behavior, is proportional to W and does not depend on U . 



4 Summary 



The low-temperature (8 = ksT) decay of the thermal diffusion in the tilted 
periodic potential (period L, amplitude Uq and external tilting force W), which 
can be described by the Fokker-Planck (FP) equation, is investigated in terms 
of the (generalized) WKB method. The nonequilibrium steady state demands 
a subtle approach to the asymptotic decay phenomenon, qualitatively different 
from the one often used in the context of the decay to the thermal equilibrium. 
However, we have confined ourselves to the manifold of k = which includes 
the steady state. This recourse makes the analysis of the decay phenomenon 
accessible. 

While eigenvalues of FP operator are complex in general, in a small tilting 
case (W < Uq) the imaginary parts of the eigenvalues are shown to be van- 
ishing. Then the Schrodinger equation associated with FP equation has the 
WKB solution satisfying PBC. Among the continuum of Bloch states, the 
zero-wavenumber eigenvalues are essential: the lowest eigevalue is responsible 
for the steady state and the decay rate is determined by the second-lowest 
eigenvalue within the zero-wavenumber manifold. The parameter dependence 
of the decay rate is explained by the tunneling process through a triple-well 

effective periodic potential. The decay rate A is given by {jf^j ^1 — Uq. 

In a large tilting case (W ^> U ), the imaginary parts of the eigenvalues 
of FP operator are crucial. We apply the complex-valued WKB method to 
the Schrodinger equation with the absorbing boundary condition, finding the 

decay phenomenon characterized by the decay rate E 1 = (j^J 8 and the 

oscillation frequency E 2 = (jf) W. The intermediate tilting case is also ex- 
plored, and our theoretical results from the WKB analysis explain W and U 
dependence of the numerical decay rate quite well for all parameter ranges 
in Fig. 3. The decay phenomena starting from k ^ distributions and the 
contribution from continuous spectra to nonexponential decays will constitute 
subjects which we intend to study in future. 
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A Perturbative analysis of the complex eigenvalue problem for 
large tilting case 

For the large tilting case (W 3> U ), the WKB quantization condition under 
the absorbing boundary condition is given by Eq.(36) or equivalently by 



We have solved Eq.(A.l) perturbatively with use of the expansion parameter 
e = ^ as follows. 

We restrict ourselves to the case n — ±1 which is observed numerically in the 
parameter range 3 < U < 20, U < W < 20 and use a branch of square-root 
with positive imaginary part in the left hand side of (A.l) (another branch 
does not satisfy (A.l)). At first, we note that (A.l) is equivalent to 




(A.l) 



o 




-L + i 



26L 



W ' 




B 



7Y 



6L 2 

^w 2 



E 2 . 



(A.2) 



We expand eigenvalue E x + iE 2 up to the second order of e as 



7T 



6L 2 
- 2 W 2 



E 1 =A 1 e + A 2 e 2 + 0(e 3 ), 



7T 



6L 2 

rW2 



E 2 = B 1 e + B 2 e 2 + 0(e 3 ). (A.3) 
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Then the real part of the equation (A. 2) is expanded in e as 



L 

/ dxll + _(_Ai + 2sin— )e 
o 

+ (^(-^ + 2 sin ^) 2 - \A 2 + I sin 2 ^ + itf 2 ) e 2 + 0( e 3 )} = L. 

(A.4) 



Noting that coefficients of each power of e should vanish and using J L c&r sin 



0, J L sin 2 ^ = ±, we obtai 



L 



am 



A 1 = 0, A 2 = ^B 2 . (A.5) 



In the same way, the imaginary part of (A. 2) is (provided A 1 = 0) 

/ «fc {f £ + i ( Sz - * sm «> + 0(£ 3)} = ^. (A . 6) 
o 



Thus one has 



46» 2 46» 

^ = 77^ 5 i = ±7T- ( A - 7 ) 



From these expressions the eigenvalue is given by Eq.(38). As shown in Fig. 3 
this estimation of E 1 is very accurate for W ^> U . 
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